Interactive ggplots

#install.packages("animint2")#from CRAN (stable).
##OR from GitHub (development).
#if(!require("remotes"))install.packages("remotes")
#remotes:install_github("tdhock/animint2")
#install.packages("remote")

library(animint2)
library(tidyverse)
library(janitor)
library(here)
library(ggmap)
library(albersusa)
library(skimr)

data(WorldBank, package="animint2")
tail (WorldBank)

##       iso2c  country year fertility.rate life.expectancy population
## 13033    ZW Zimbabwe 2007          3.491        45.79707   12740160
## 13034    ZW Zimbabwe 2008          3.428        47.07061   12784041
## 13035    ZW Zimbabwe 2009          3.360        48.45049   12888918
## 13036    ZW Zimbabwe 2010          3.290        49.86088   13076978
## 13037    ZW Zimbabwe 2011          3.219        51.23644   13358738
## 13038    ZW Zimbabwe 2012             NA              NA   13724317
##       GDP.per.capita.Current.USD 15.to.25.yr.female.literacy iso3c
## 13033                   415.3755                          NA   ZWE
## 13034                   345.4074                          NA   ZWE
## 13035                   475.8538                          NA   ZWE
## 13036                   568.4275                    99.55316   ZWE
## 13037                   722.8377                          NA   ZWE
## 13038                   787.9382                          NA   ZWE
##                                       region capital longitude latitude
## 13033 Sub-Saharan Africa (all income levels)  Harare   31.0672 -17.8312
## 13034 Sub-Saharan Africa (all income levels)  Harare   31.0672 -17.8312
## 13035 Sub-Saharan Africa (all income levels)  Harare   31.0672 -17.8312
## 13036 Sub-Saharan Africa (all income levels)  Harare   31.0672 -17.8312
## 13037 Sub-Saharan Africa (all income levels)  Harare   31.0672 -17.8312
## 13038 Sub-Saharan Africa (all income levels)  Harare   31.0672 -17.8312
##           income lending
## 13033 Low income   Blend
## 13034 Low income   Blend
## 13035 Low income   Blend
## 13036 Low income   Blend
## 13037 Low income   Blend
## 13038 Low income   Blend

WorldBank1975 <- subset(WorldBank, year==1975)
head(WorldBank1975)

##     iso2c              country year fertility.rate life.expectancy population
## 281    AD              Andorra 1975             NA              NA      30706
## 334    AE United Arab Emirates 1975          6.009        66.18539     532742
## 387    AF          Afghanistan 1975          7.692        37.25712   12551790
## 440    AG  Antigua and Barbuda 1975             NA              NA      69253
## 493    AL              Albania 1975          4.417        68.32583    2426592
## 546    AM              Armenia 1975          2.745        70.52751    2825650
##     GDP.per.capita.Current.USD 15.to.25.yr.female.literacy iso3c
## 281                  7168.3987                          NA   AND
## 334                 27631.8985                    56.27697   ARE
## 387                   188.5521                          NA   AFG
## 440                         NA                          NA   ATG
## 493                         NA                          NA   ALB
## 546                         NA                          NA   ARM
##                                             region          capital longitude
## 281      Europe & Central Asia (all income levels) Andorra la Vella    1.5218
## 334 Middle East & North Africa (all income levels)        Abu Dhabi   54.3705
## 387                                     South Asia            Kabul   69.1761
## 440  Latin America & Caribbean (all income levels)     Saint John's  -61.8456
## 493      Europe & Central Asia (all income levels)           Tirane   19.8172
## 546      Europe & Central Asia (all income levels)          Yerevan    44.509
##     latitude               income        lending
## 281  42.5075 High income: nonOECD Not classified
## 334  24.4764 High income: nonOECD Not classified
## 387  34.5228           Low income            IDA
## 440  17.1175  Upper middle income           IBRD
## 493  41.3317  Upper middle income           IBRD
## 546  40.1596  Lower middle income          Blend

scatter1 <- ggplot()+
  geom_point(
    mapping=aes(x=life.expectancy, y=fertility.rate, color=region),
    data=WorldBank1975)
scatter1

## Warning: Removed 27 rows containing missing values (geom_point).

#animint(scatter1)

WorldBankBefore1975 <- subset(WorldBank, 1970 <= year & year <= 1975)
two.layers <- scatter1+
  geom_path(aes(x=life.expectancy, y=fertility.rate, color=region,
                group=country),
            data=WorldBankBefore1975)

two.layers

## Warning: Removed 27 rows containing missing values (geom_point).

## Warning: Removed 162 rows containing missing values (geom_path).

#(viz.two.layers <- animint(two.layers)) Generates an error: 'unitType' is not an exported object from 'namespace:grid'

three.layers <- two.layers+
  geom_text(aes(x=life.expectancy, y=fertility.rate, color=region,
                label=country),
            data=WorldBank1975)

three.layers

## Warning: Removed 27 rows containing missing values (geom_point).

## Warning: Removed 162 rows containing missing values (geom_path).

## Warning: Removed 27 rows containing missing values (geom_text).

two.layersstorage <- two.layers
two.layers$timeSeries <- ggplot()+
  geom_line(aes(x=year, y=fertility.rate, color=region, group=country),
            data=WorldBank)

two.layersstorage

## Warning: Removed 27 rows containing missing values (geom_point).

## Warning: Removed 162 rows containing missing values (geom_path).

two.layers$timeSeries

## Warning: Removed 1183 rows containing missing values (geom_path).

add.x.var <- function(df, x.var){
  data.frame(df, x.var=factor(x.var, c("life expectancy", "year")))
}

  scatterfirst=ggplot()+
    geom_point(aes(x=life.expectancy, y=fertility.rate, color=region),
               data=add.x.var(WorldBank1975, "life expectancy"))+
    geom_path(aes(x=life.expectancy, y=fertility.rate, color=region,
                  group=country),
              data=add.x.var(WorldBankBefore1975, "life expectancy"))+
    geom_line(aes(x=year, y=fertility.rate, color=region, group=country),
              data=add.x.var(WorldBank, "year"))+theme(axis.text.x = element_text(angle = 50, hjust = 1))+
    xlab("")+
    facet_grid(. ~ x.var, scales="free")+
    theme_bw()+
    theme(panel.margin=grid::unit(0, "lines"))

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

  scatterfirst

## Warning: Removed 27 rows containing missing values (geom_point).

## Warning: Removed 162 rows containing missing values (geom_path).

## Warning: Removed 1183 rows containing missing values (geom_path).

show.point.list <- list()
show.path.list <- list()
for(show.year in c(1975, 1985, 1995)){
  show.point.list[[paste(show.year)]] <- data.frame(
    show.year, subset(WorldBank, year==show.year))
  show.path.list[[paste(show.year)]] <- data.frame(
    show.year, subset(WorldBank, show.year - 5 <= year & year <= show.year))
}
show.point <- do.call(rbind, show.point.list)
show.path <- do.call(rbind, show.path.list)

threepanelsvisual <- list(
  scatter=ggplot()+
    geom_point(aes(x=life.expectancy, y=fertility.rate, color=region),
               data=show.point)+
    geom_path(aes(x=life.expectancy, y=fertility.rate, color=region,
                  group=country),
              data=show.path)+
    facet_grid(. ~ show.year)+
    theme_bw()+
    theme(panel.margin=grid::unit(0, "lines"))
  )

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

threepanelsvisual

## $scatter

## Warning: Removed 76 rows containing missing values (geom_point).

## Warning: Removed 422 rows containing missing values (geom_path).

Exercises

Visualize a different type of data not present in the animint2 Manual.

Data Source: https://rcweb.dartmouth.edu/~f002d69/workshops/index_rspatial.html which represents sightings of grizzly bears (Ursos arctos) in Alaska.

Other data set used is from Brandell, Ellen E (2021), Serological dataset and R code for: Patterns and processes of pathogen exposure in gray wolves across North America, Dryad, (https://doi.org/10.5061/dryad.5hqbzkh51).

Load the grizzly data and evaluate its structure.

grizzly <- read_csv(here("data", "bear-sightings.csv")) 

## Parsed with column specification:
## cols(
##   bear.id = col_double(),
##   longitude = col_double(),
##   latitude = col_double()
## )

grizzly <- janitor::clean_names(grizzly)
grizzly

## # A tibble: 494 x 3
##    bear_id longitude latitude
##      <dbl>     <dbl>    <dbl>
##  1       7     -149.     62.7
##  2      57     -153.     58.4
##  3      69     -145.     62.4
##  4      75     -153.     59.9
##  5     104     -143.     61.1
##  6     108     -150.     62.9
##  7     115     -152.     68.0
##  8     116     -147.     62.6
##  9     125     -157.     60.2
## 10     135     -156.     58.9
## # … with 484 more rows

glimpse(grizzly)

## Rows: 494
## Columns: 3
## $ bear_id   <dbl> 7, 57, 69, 75, 104, 108, 115, 116, 125, 135, 137, 162, 185,…
## $ longitude <dbl> -148.9560, -152.6228, -144.9374, -152.8485, -143.2948, -149…
## $ latitude  <dbl> 62.65822, 58.35064, 62.38227, 59.90122, 61.07311, 62.91605,…

summary(grizzly)

##     bear_id       longitude         latitude    
##  Min.   :   7   Min.   :-166.2   Min.   :55.02  
##  1st Qu.:2569   1st Qu.:-154.2   1st Qu.:58.13  
##  Median :4822   Median :-151.0   Median :60.97  
##  Mean   :4935   Mean   :-149.1   Mean   :61.41  
##  3rd Qu.:7387   3rd Qu.:-145.6   3rd Qu.:64.13  
##  Max.   :9996   Max.   :-131.3   Max.   :70.37

str(grizzly)

## tibble [494 × 3] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ bear_id  : num [1:494] 7 57 69 75 104 108 115 116 125 135 ...
##  $ longitude: num [1:494] -149 -153 -145 -153 -143 ...
##  $ latitude : num [1:494] 62.7 58.4 62.4 59.9 61.1 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   bear.id = col_double(),
##   ..   longitude = col_double(),
##   ..   latitude = col_double()
##   .. )

skim(grizzly)

Data summary

Name	grizzly
Number of rows	494
Number of columns	3
_______________________
Column type frequency:
numeric	3
________________________
Group variables	None

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
bear_id	0	1	4934.51	2877.87	7.00	2568.75	4822.50	7386.75	9996.00	▇▇▇▆▇
longitude	0	1	-149.08	7.65	-166.23	-154.25	-150.96	-145.59	-131.26	▁▇▇▂▃
latitude	0	1	61.41	3.83	55.02	58.13	60.97	64.13	70.37	▇▇▇▃▃

Use the range of the latitude and longitude to build an appropriate bounding box for a map.

latitude_me <- c(55.02,70.37)
longitude_me <- c(-166.2, -131.3)
bear_box <- make_bbox(longitude_me, latitude_me, f=0.05)

Load a map from stamen in a terrain style projection and display the map.

map_projection <- get_map(bear_box, maptype = "terrain", source = "stamen")

## Source : http://tile.stamen.com/terrain/5/1/6.png

## Source : http://tile.stamen.com/terrain/5/2/6.png

## Source : http://tile.stamen.com/terrain/5/3/6.png

## Source : http://tile.stamen.com/terrain/5/4/6.png

## Source : http://tile.stamen.com/terrain/5/1/7.png

## Source : http://tile.stamen.com/terrain/5/2/7.png

## Source : http://tile.stamen.com/terrain/5/3/7.png

## Source : http://tile.stamen.com/terrain/5/4/7.png

## Source : http://tile.stamen.com/terrain/5/1/8.png

## Source : http://tile.stamen.com/terrain/5/2/8.png

## Source : http://tile.stamen.com/terrain/5/3/8.png

## Source : http://tile.stamen.com/terrain/5/4/8.png

## Source : http://tile.stamen.com/terrain/5/1/9.png

## Source : http://tile.stamen.com/terrain/5/2/9.png

## Source : http://tile.stamen.com/terrain/5/3/9.png

## Source : http://tile.stamen.com/terrain/5/4/9.png

## Source : http://tile.stamen.com/terrain/5/1/10.png

## Source : http://tile.stamen.com/terrain/5/2/10.png

## Source : http://tile.stamen.com/terrain/5/3/10.png

## Source : http://tile.stamen.com/terrain/5/4/10.png

ggmap(map_projection)

Build a final map that overlays the recorded observations of grizzly bears in Alaska.

ggmap(map_projection) + geom_point(data = grizzly, aes(longitude,latitude), size=1, color = "blue") + labs(x = "Longitude", y= "Latitude", title = "Recorded Observations of Grizzly Bears in Alaska")

Exercise 2: Wolves data. Load the data and evaluate its structure.

wolves <- readr::read_csv("data/wolves_data/wolves_dataset.csv")

## Parsed with column specification:
## cols(
##   .default = col_double(),
##   pop = col_character(),
##   age.cat = col_character(),
##   sex = col_character(),
##   color = col_character()
## )

## See spec(...) for full column specifications.

str(wolves)

## tibble [1,986 × 23] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ pop               : chr [1:1986] "AK.PEN" "AK.PEN" "AK.PEN" "AK.PEN" ...
##  $ year              : num [1:1986] 2006 2006 2006 2006 2006 ...
##  $ age.cat           : chr [1:1986] "S" "S" "A" "S" ...
##  $ sex               : chr [1:1986] "F" "M" "F" "M" ...
##  $ color             : chr [1:1986] "G" "G" "G" "B" ...
##  $ lat               : num [1:1986] 57 57 57 57 57 ...
##  $ long              : num [1:1986] -158 -158 -158 -158 -158 ...
##  $ habitat           : num [1:1986] 254 254 254 254 254 ...
##  $ human             : num [1:1986] 10.4 10.4 10.4 10.4 10.4 ...
##  $ pop.density       : num [1:1986] 8 8 8 8 8 8 8 8 8 8 ...
##  $ pack.size         : num [1:1986] 8.78 8.78 8.78 8.78 8.78 8.78 8.78 8.78 8.78 8.78 ...
##  $ standard.habitat  : num [1:1986] -1.63 -1.63 -1.63 -1.63 -1.63 ...
##  $ standard.human    : num [1:1986] -0.978 -0.978 -0.978 -0.978 -0.978 ...
##  $ standard.pop      : num [1:1986] -0.683 -0.683 -0.683 -0.683 -0.683 ...
##  $ standard.packsize : num [1:1986] 1.32 1.32 1.32 1.32 1.32 ...
##  $ standard.latitude : num [1:1986] 0.721 0.721 0.721 0.721 0.721 ...
##  $ standard.longitude: num [1:1986] -2.14 -2.14 -2.14 -2.14 -2.14 ...
##  $ cav.binary        : num [1:1986] 1 1 1 0 1 1 1 1 1 1 ...
##  $ cdv.binary        : num [1:1986] 0 0 0 0 0 1 0 0 0 0 ...
##  $ cpv.binary        : num [1:1986] 0 0 1 1 0 1 0 0 0 0 ...
##  $ chv.binary        : num [1:1986] 1 1 1 1 1 1 1 1 1 0 ...
##  $ neo.binary        : num [1:1986] NA NA NA 0 0 NA NA 1 0 1 ...
##  $ toxo.binary       : num [1:1986] NA NA NA 1 0 NA NA 1 0 0 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   pop = col_character(),
##   ..   year = col_double(),
##   ..   age.cat = col_character(),
##   ..   sex = col_character(),
##   ..   color = col_character(),
##   ..   lat = col_double(),
##   ..   long = col_double(),
##   ..   habitat = col_double(),
##   ..   human = col_double(),
##   ..   pop.density = col_double(),
##   ..   pack.size = col_double(),
##   ..   standard.habitat = col_double(),
##   ..   standard.human = col_double(),
##   ..   standard.pop = col_double(),
##   ..   standard.packsize = col_double(),
##   ..   standard.latitude = col_double(),
##   ..   standard.longitude = col_double(),
##   ..   cav.binary = col_double(),
##   ..   cdv.binary = col_double(),
##   ..   cpv.binary = col_double(),
##   ..   chv.binary = col_double(),
##   ..   neo.binary = col_double(),
##   ..   toxo.binary = col_double()
##   .. )

skim(wolves)

Data summary

Name	wolves
Number of rows	1986
Number of columns	23
_______________________
Column type frequency:
character	4
numeric	19
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
pop	0	1.00	2	7	0	17	0
age.cat	102	0.95	1	1	0	3	0
sex	19	0.99	1	1	0	2	0
color	616	0.69	1	1	0	2	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
year	0	1.00	2010.28	5.91	1992.00	2006.00	2011.00	2016.00	2019.00	▁▂▆▆▇
lat	0	1.00	50.43	9.16	33.89	44.60	46.83	57.89	80.50	▁▇▂▂▁
long	0	1.00	-116.86	19.11	-157.84	-123.73	-110.99	-110.55	-82.42	▃▁▂▇▂
habitat	0	1.00	12797.42	7330.64	254.08	10375.16	11211.30	11860.81	34676.55	▂▇▂▁▁
human	0	1.00	2335.38	1912.21	0.02	80.60	2787.67	3973.47	6228.64	▇▁▅▅▁
pop.density	0	1.00	14.91	9.50	3.74	7.40	11.63	25.32	33.96	▇▅▁▂▅
pack.size	0	1.00	6.47	1.67	3.55	5.62	6.37	8.25	9.56	▅▇▅▆▃
standard.habitat	0	1.00	0.01	0.96	-1.63	-0.31	-0.20	-0.11	2.88	▂▇▂▁▁
standard.human	0	1.00	0.15	0.92	-0.98	-0.94	0.36	0.94	2.03	▇▁▅▅▁
standard.pop	0	1.00	0.05	1.01	-1.13	-0.75	-0.30	1.15	2.07	▇▅▁▂▅
standard.packsize	0	1.00	-0.04	0.98	-1.76	-0.54	-0.10	1.00	1.77	▅▇▅▆▃
standard.latitude	0	1.00	0.00	1.00	-1.81	-0.64	-0.39	0.81	3.28	▁▇▂▂▁
standard.longitude	0	1.00	0.00	1.00	-2.14	-0.36	0.31	0.33	1.80	▃▁▂▇▂
cav.binary	321	0.84	0.85	0.35	0.00	1.00	1.00	1.00	1.00	▂▁▁▁▇
cdv.binary	21	0.99	0.22	0.42	0.00	0.00	0.00	0.00	1.00	▇▁▁▁▂
cpv.binary	7	1.00	0.79	0.40	0.00	1.00	1.00	1.00	1.00	▂▁▁▁▇
chv.binary	548	0.72	0.80	0.40	0.00	1.00	1.00	1.00	1.00	▂▁▁▁▇
neo.binary	538	0.73	0.28	0.45	0.00	0.00	0.00	1.00	1.00	▇▁▁▁▃
toxo.binary	827	0.58	0.48	0.50	0.00	0.00	0.00	1.00	1.00	▇▁▁▁▇

str(wolves)

## tibble [1,986 × 23] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ pop               : chr [1:1986] "AK.PEN" "AK.PEN" "AK.PEN" "AK.PEN" ...
##  $ year              : num [1:1986] 2006 2006 2006 2006 2006 ...
##  $ age.cat           : chr [1:1986] "S" "S" "A" "S" ...
##  $ sex               : chr [1:1986] "F" "M" "F" "M" ...
##  $ color             : chr [1:1986] "G" "G" "G" "B" ...
##  $ lat               : num [1:1986] 57 57 57 57 57 ...
##  $ long              : num [1:1986] -158 -158 -158 -158 -158 ...
##  $ habitat           : num [1:1986] 254 254 254 254 254 ...
##  $ human             : num [1:1986] 10.4 10.4 10.4 10.4 10.4 ...
##  $ pop.density       : num [1:1986] 8 8 8 8 8 8 8 8 8 8 ...
##  $ pack.size         : num [1:1986] 8.78 8.78 8.78 8.78 8.78 8.78 8.78 8.78 8.78 8.78 ...
##  $ standard.habitat  : num [1:1986] -1.63 -1.63 -1.63 -1.63 -1.63 ...
##  $ standard.human    : num [1:1986] -0.978 -0.978 -0.978 -0.978 -0.978 ...
##  $ standard.pop      : num [1:1986] -0.683 -0.683 -0.683 -0.683 -0.683 ...
##  $ standard.packsize : num [1:1986] 1.32 1.32 1.32 1.32 1.32 ...
##  $ standard.latitude : num [1:1986] 0.721 0.721 0.721 0.721 0.721 ...
##  $ standard.longitude: num [1:1986] -2.14 -2.14 -2.14 -2.14 -2.14 ...
##  $ cav.binary        : num [1:1986] 1 1 1 0 1 1 1 1 1 1 ...
##  $ cdv.binary        : num [1:1986] 0 0 0 0 0 1 0 0 0 0 ...
##  $ cpv.binary        : num [1:1986] 0 0 1 1 0 1 0 0 0 0 ...
##  $ chv.binary        : num [1:1986] 1 1 1 1 1 1 1 1 1 0 ...
##  $ neo.binary        : num [1:1986] NA NA NA 0 0 NA NA 1 0 1 ...
##  $ toxo.binary       : num [1:1986] NA NA NA 1 0 NA NA 1 0 0 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   pop = col_character(),
##   ..   year = col_double(),
##   ..   age.cat = col_character(),
##   ..   sex = col_character(),
##   ..   color = col_character(),
##   ..   lat = col_double(),
##   ..   long = col_double(),
##   ..   habitat = col_double(),
##   ..   human = col_double(),
##   ..   pop.density = col_double(),
##   ..   pack.size = col_double(),
##   ..   standard.habitat = col_double(),
##   ..   standard.human = col_double(),
##   ..   standard.pop = col_double(),
##   ..   standard.packsize = col_double(),
##   ..   standard.latitude = col_double(),
##   ..   standard.longitude = col_double(),
##   ..   cav.binary = col_double(),
##   ..   cdv.binary = col_double(),
##   ..   cpv.binary = col_double(),
##   ..   chv.binary = col_double(),
##   ..   neo.binary = col_double(),
##   ..   toxo.binary = col_double()
##   .. )

Evaluating how many distinct wolf populations are included in this study:

wolves$pop <- as.factor(wolves$pop)
n_distinct(wolves$pop)

## [1] 17

wolves_US <- wolves %>%
  filter(lat<=48)
wolves_US

## # A tibble: 1,169 x 23
##    pop    year age.cat sex   color   lat  long habitat human pop.density
##    <fct> <dbl> <chr>   <chr> <chr> <dbl> <dbl>   <dbl> <dbl>       <dbl>
##  1 GTNP   2012 P       M     G      43.8 -111.  10375. 3924.        34.0
##  2 GTNP   2012 P       F     G      43.8 -111.  10375. 3924.        34.0
##  3 GTNP   2012 P       F     G      43.8 -111.  10375. 3924.        34.0
##  4 GTNP   2012 P       M     B      43.8 -111.  10375. 3924.        34.0
##  5 GTNP   2013 A       F     G      43.8 -111.  10375. 3924.        34.0
##  6 GTNP   2013 A       M     G      43.8 -111.  10375. 3924.        34.0
##  7 GTNP   2013 P       M     G      43.8 -111.  10375. 3924.        34.0
##  8 GTNP   2013 P       M     G      43.8 -111.  10375. 3924.        34.0
##  9 GTNP   2013 P       M     G      43.8 -111.  10375. 3924.        34.0
## 10 GTNP   2013 P       F     G      43.8 -111.  10375. 3924.        34.0
## # … with 1,159 more rows, and 13 more variables: pack.size <dbl>,
## #   standard.habitat <dbl>, standard.human <dbl>, standard.pop <dbl>,
## #   standard.packsize <dbl>, standard.latitude <dbl>, standard.longitude <dbl>,
## #   cav.binary <dbl>, cdv.binary <dbl>, cpv.binary <dbl>, chv.binary <dbl>,
## #   neo.binary <dbl>, toxo.binary <dbl>

Use the albersusa package to make a base map of the lower 48 US states.

us_lower <- usa_sf()

c_comp <- counties_sf()

ggplot() +
  geom_sf(data=us_lower, size = 0.25) + theme_linedraw() + labs(title = "US State Boundaries")

Use the relimited data to plot the distribution of wolf populations in the lower 48 US states.

ggplot() + 
  geom_sf(data = us_lower, size = 0.125) + 
  theme_linedraw()+
  geom_point(data = wolves_US, aes(long, lat, color=pop), shape=2, size=4)+
  labs(title = "Distribution of Wolf Species",
       x = "Longitude",
       y = "Latitude")+
  theme_minimal()

Average pack size for the wolves in this study by region:

wolves %>%
  group_by(pop) %>%
  summarize(avg_wolfpacksize = mean(pack.size), .groups='keep') %>%
  arrange(desc(avg_wolfpacksize))

## # A tibble: 17 x 2
## # Groups:   pop [17]
##    pop     avg_wolfpacksize
##    <fct>              <dbl>
##  1 BAN.JAS             9.56
##  2 ELLES               9.19
##  3 AK.PEN              8.78
##  4 YNP                 8.25
##  5 GTNP                8.1 
##  6 MI                  7.12
##  7 DENALI              6.45
##  8 YUCH                6.37
##  9 INT.AK              6.24
## 10 BC                  5.88
## 11 MT                  5.62
## 12 SE.AK               5   
## 13 SNF                 4.81
## 14 ONT                 4.37
## 15 MEXICAN             4.04
## 16 N.NWT               4   
## 17 SS.NWT              3.55

Make a new map that shows the distribution of wolves in the lower 48 US states but which has the size of location markers adjusted by pack size.

wolfplot <- ggplot() +
  geom_sf(data=us_lower, size=0.2) + 
  geom_point(data=wolves_US, aes(long,lat,size=pack.size,color=pack.size), shape=5) + theme_minimal()+
  labs(title = "Distribution of Wolves in Lower 48 US states by Pack Size", x = "Longitude", y="Latitude")

wolfplot